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SUMMARY 

For convection-dominated flows, classical second-order methods are notoriously 
oscillatory and often unstable. For this reason, many computational-fluid-dynam- 
icists have adopted various forms of (inherently stable) first-order upwinding over 
the past few decades. Although it is now well known that first-order convection 
schemes suffer from serious inaccuracies attributable to artificial viscosity or 
numerical diffusion under high-convection conditions, these methods continue to 
enjoy widespread popularity for numerical heat-transfer calculations, apparently 
due to a perceived lack of viable high-accuracy alternatives. But alternatives are 
available. For example, nonoscillatory methods used in gasdynamics, including 
currently popular "TVD” schemes, can be easily adapted to multidimensional incom- 
pressible flow and convective transport. This, in itself, would be a major advance for 
numerical convective heat transfer, for example. But, as this report shows, second- 
order TVD schemes form only a small, overly restrictive, subclass of a much more 
universal, and extremely simple, nonoscillatory flux-limiting strategy which can be 
applied to convection schemes of arbitrarily high-order accuracy, while requiring 
only a simple tridiagonal ADI line-solver, as used in the majority of general-purpose 
iterative codes for incompressible flow and numerical heat transfer. The new uni- 
versal limiter and associated solution procedures form the so-called ULTRA-SHARP 
alternative for high-resolution nonoscillatory multidimensional steady-state high- 
speed convective modelling. 


INTRODUCTION 

For many years the state of the art in high-speed convective modelling, espe- 
cially in the field of numerical heat and mass transfer, has been dominated by first- 
order upwinding, often in the guise of the "Hybrid” scheme of Spalding [1] or, more 
recently, the related "Power-Law Differencing Scheme” (PLDS) of Patankar [2]. 
This situation has clearly evolved from an attempt to remedy the infamous problems 
of unphysical oscillations and instabilities associated with "classical” central-differ- 
ence methods under high-convection conditions, using practical grids. The Hybrid 
scheme avoids oscillatory behaviour by switching from second-order central to first- 
order upwinding for convection (and omitting modelled physical diffusion) wherever 
the local component grid Peclet (or Reynolds) number exceeds a value of 2. PLDS 
and the exponential differencing scheme (EDS) on which it is based [3] involve a 
more subtle blending strategy, but both are also equivalent to first-order upwinding 
for convection (with physical diffusion totally suppressed) for component grid Peclet 
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numbers greater than about 6. For most cases of practical interest, especially in 
three dimensions where very fine meshes are out of the question, the actual (laminar 
or turbulent) component grid Peclet or Reynolds numbers are likely to be orders of 
magnitude larger than 2 or 6 throughout the bulk of the flow domain. This means 
that Hybrid and PLDS are operating as first-order upwinding almost everywhere in 
the flow-field except at a very small fraction of grid points where the flow direction 
happens to be parallel to grid lines or near boundaries and stagnation regions where 
the convecting velocities are small. Thus, instead of solving a high-convection 
problem, these methods simulate an unphysical (and anisotropic) Zoio-convection 
problem in which the effective local component grid Peclet or Reynolds number can 
never be greater than 2. 

As has been well established theoretically for over a decade [4-6], first-order 
upwinding for convection is a highly inaccurate method because of the effects of 
global artificial numerical diffusion (or artificial viscosity). Thus, in many cases of 
practical interest, Hybrid or PLDS may give rather low-accuracy results; this has 
been directly confirmed in several comparative benchmark studies [7-10], Hybrid 
and PLDS, as used in various evolutions of the well-known TEACH code [11], have 
been most successful for flows in which the convective terms are not important; i.e., 
for essentially potential flow governed by the kinematic constraint of the continuity 
equation (V • v = 0), modified by boundary regions which are usually handled by 
built-in wall-functions , and for turbulence modelling when turbulent-transport 
equations are dominated by source-sink terms [12]. However, for simulating recircu- 
lating flows or strong shear flows or for predicting scalar transport, these methods 
must be considered to be highly unreliable because of the inherent global artificial 
viscosity or diffusivity. 

Even though these defects are well documented and generally acknowledged, 
the Hybrid and PLDS methods continue to enjoy widespread popularity, apparently 
because available alternatives claiming better accuracy (such as higher-order 
upwinding) are either not well known or have been discarded because of perceived 
incompatibilities with TEACH-like solution algorithms [10]. When carefully 
applied to three-dimensional laminar recirculating flows however, higher-order 
upwinding methods such a QUICK [13] have been extremely successful [14], show- 
ing fine-scale details on relatively coarse grids that are simply wiped out by artificial 
numerical viscosity in the corresponding Hybrid calculation [15]. In many convec- 
tion-dominated flow problems, the multidimensional third-order QUICK scheme has 
been shown to be stable, economical, and highly accurate; because of the wider com- 
putational stencil, the solution algorithm is most naturally cast in penfadiagonal 
form [16], but can also be used with conventional tridiagonal solvers provided that 
outlying and other compensating terms [17] are carefully transferred to the source 
term. 


However, in some cases (such as turbulent flows, for example) there may be 
problems with higher-order upwinding because, in regions involving sudden jumps 
in value, there is a tendency for the simulation to overshoot or undershoot the correct 
transition value. In many cases, this presents no problem other than slight inaccu- 
racy; however, in codes where turbulent transport variables (such as eddy viscosity, 
for example) are computed as part of the solution procedure, it is possible for a 
modelled transport coefficient to undershoot to a point where it takes on locally 
negative values, thus resulting in violent nonlinear instability. Similar defects 
occur in various forms of skew upwinding [18,8]. There have been recent attempts to 
suppress the undershoot problem by blending higher-order (or skew) upwinding 
methods with standard (component-wise) first-order upwinding, using sophisticated 
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blending strategies [19,20], Although much more accurate than Hybrid or PLDS, 
these methods still generate an unnecessary degree of numerical smearing when 
attempting to simulate near-discontinuities. 


The current situation in numerical gasdynamics is considerably better. The 
most popular schemes are based on second-order central differencing for all spatia 
derivatives, as used in the canonical second-order Lax-Wendroff method 121], lor 
example, together with explicitly added global fourth-order artificial dissipation (a 
spatial fourth-difference term) to enhance convergence in smooth regions and an 
explicit locally varying nonlinear solution-dependent artificial difrusivity or viscos- 
ity (a spatial second-difference) to inhibit oscillations in regions of sharply varying 
gradient [22,23]. This explicit damping of a potentially oscillatory scheme is usually 
quite successful in giving second-order accuracy without excessive oscillation or 
gross artificial smearing. In general terms, Hybrid and PLDS could be considered to 
follow a similar philosophy: i.e., using second-order central differencing with locally 
varying solution-dependent artificial diffusivity or viscosity added explicitly . But a 
major difference lies in the strategy adopted for the evaluation of the magnitude or 
the artificial terms. With Hybrid and PLDS it is a function of the local component 
grid Peclet or Reynolds number, so that artificial diffusivity is always large 
wherever physical convection dominates physical diffusion at the grid-cell level, 
thus guaranteeing that modelled convection never dominates numerical diffusion. 
By contrast, the magnitude of the gasdynamics codes’ artificial diffusivity depends in 
a nonlinear way on local behaviour of the convected field variables: very low in 
smooth regions, but rapidly increasing in regions of strongly varying gradients (high 
curvature), irrespective of the local Reynolds number, which could be infinite (as in 
Euler-equation simulations). Thus, in these schemes, artificial viscosity is explicitly 
added (automatically) only where "needed” (to suppress potential oscillations near 
shock waves or contact discontinuities, for example). This technique^has been devel- 
oped into a sophisticated art with so-called "shock-capturing” and "total-variation- 
diminishing” (TVD) schemes [24,25], which may even involve local negative arti- 
ficial viscosity, resulting in "artificial compression” to aid in numerically steepening 
discontinuities without overshoot. The popular "Superbee” scheme is of this latter 
type [26]. 


Apparently there has never been any attempt to explicitly adapt the highly 
successful second-order gasdynamic codes or the so-called "high-resolution (but 
formally still second-order) TVD schemes to incompressible flows and convective 
heat and mass transfer, even though (as will be shown) this is extremely simple to 
do, using either explicit time-marching or an iterative pentadiagonal matrix 
algorithm based on a variable-curvature- factor technique [27,28]. But it turns out 
that the commonly used TVD flux limiters form an overly restrictive subclass ol a 
much more universal technique guaranteeing tight resolution of discontinuities 
without overshoots or undershoots; and this "universal limiter can be applied to 
methods of arbitrarily high-order accuracy, giving very sharp resolution without 
introducing artificial numerical compression (and concomitant distortion and 
clipping of smooth profiles [29]). The strategy developed here is based on a simple 
high-accuracy resolution program which uses third-order upwinding in smooth 
regions and adaptively increases the accuracy (by locally expanding the compu- 
tational stencil automatically) in high-curvature regions, while applying the 
universal limiter to eliminate any possibility of overshoot or oscillation. Thus, the 
Universal Limiter for Tight Resolution and Accuracy in combination with this 
Simple High-Accuracy Resolution Program results in the "ULTRA-SHARP alter- 
native for high-speed nonoscillatory steady-state convective modelling. 
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mlrlr t damor ? str ated on a number of simple two-dimensional steady-state bench- 
mark test problems for scalar transport in prescribed velocity fields. 


SHORT COMINGS OF HYBRID AND PLDS 

Effective Grid Peclet Number 

for a sca"a S r d r **“ m ° de ' convection ' dilT nsion transport equation in one dimension 


AJ> 

dt 


<M> d% 
u— + D — 


dx 


dx * 


(1) 


for constant u and D. The forward-time, central 
equation takes the form [30] 


space (FTCS) discretization of this 


4>" + 1 - 4," 

At 


= — u 


( *; +1 - ( 4 > r +1 - 2 < +** ,) 

+ D ? _ 


2 Ax 


Ax* 


(2) 


using standard index notation on a uniform grid. This can be rearranged 


as 


<t> 


n+ 1 


= *■ - ? . - *1 ,> + r «-** , - 2*" + c ,i 


introducing the Courant number 


(3) 


c — u At I Ax 

and the cell Peclet number (or Reynolds number 
vorticity transport). 


(4) 


in the case of momentum or 


P A = \u\Ax/D 


(5) 


Equation (1) can also be discretized using first-order upwindine for the convec 
tion term and simply ignoring the physical diffusion term, giving (fo? u > 0) 


<t> 


n + 1 


4>; - d 4>; - <{>;_,) 


(6) 
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with a similar form for u < 0. Equation (6) can be rewritten in the form of Equation 
(3) as 


*"*' = *r- f <+“*, - **.i» + 7 - 2 i>r + ♦?-,> 

P A 


where the effective grid Peclet number is, in this case 


P A - 2 


(7) 


(8) 


Thus, the Hybrid scheme - defined by using Equation (3) for P. ^ 2 and Equation 
(6) for P a > 2 - can be written in the form of Equation (7) with 


P. = P. for P A <; 2 
A A A 


(9) 


and 


Pi - 2 for P A > 2 (10) 

Now consider an exact solution of the steady-state form of Equation (1) with a 
downstream boundary condition <J>(L) = 1, and with <|> = 0 far upstream. This can 
be written 


4 >U) = exp[ — (L— *)/A] for i S t 

where the length-constant A is given by 


( 11 ) 


X = D/u ( 12 ) 

The corresponding central-difference equation is 

(2 — P A )< ^ > i + l — 4 ( 2 + P^) 4 >j_i = 0 ( 13 ) 

and the analytical solution of this difference equation for the above boundary condi- 
tions is 



where k is an integer defined, for discrete values of x, by 


( 14 ) 


k = (L — x)/ Ax for x < L ( 15 ) 

Note that the exact solution of the differential equation, Equation (11), can be 
written, for discrete node values, as 
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(16) 


*k 



Thus, for P A < 2, the central-difference solution is like a solution to the differential 
equation, but at a larger grid Peclet number. And for P. > 2, Equation (14) shows 
that the discrete solution becomes oscillatory, as is well known [30]. These features 
are demonstrated in Figure 1, which shows the factors in Equations (14) and (16) as a 
function of P A . 


Now consider Equation (13) rewritten in terms of an effective grid Peclet 
number, P A * : 


(2-P*)4» i + 1 -4<t>. + (2 + P*)<j > j _ 1 = ° (17) 

The solution of this equation is of the same form as Equation (14), with P A replaced 
by P A *. But this can be forced to match the true solution. Equation (16), by setting 


2 - P. 




2 + P, 


(18) 


giving, on rearrangement, an explicit formula for the effective grid Peclet number in 
terms of the physical value 


P A = 2 lanh(P A /2) (19) 

This is the basis of the exponential difference scheme [3] or so-called "optimal” 
upwinding [31]. Using P. -dependent weighting factors for convection and diffusion, 
Raithby and Schneider [3] develop algebraic approximations (to avoid expensive 
exponential evaluations) which can easily be shown to be equivalent to an approxi- 
mation to Equation (19) given by 


P A (l + 0.005 P A ) 



2(5 + P A ) P (1 + 0.05 P A ) (20) 

AAA 

for all (positive) values of P A . Similarly, Patankar’s power-law scheme [2] can be 
interpreted as 


P A (R&S) = 


2P 

P*(PLDS) = 0<P <10 

P A + 2(1 - 0.1PJ 5 

A A 

= 2 P A > 10 

A 

In the same spirit, Spalding’s Hybrid scheme [1], given by Equations (9) and (10), 
could be interpreted as a (much less accurate) piece- wise linear approximation of the 
hyperbolic-tangent function. Equations (9) and (10), (19), (20), and (21) are por- 
trayed in Figure 2. As defined, Hybrid is equivalent to first-order upwinding for 
convection with physical diffusion omitted (i.e., P A * = 2) for P A > 2. But also note 
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that the other schemes are effectively equivalent to this for P. greater than about 6. 
This is an important point which is sometimes lost in the details of the algebraic or 
power-law approximations. In other words, although (as seen from Figure 2) Equa- 
tions (20) and (21) are both very good engineering approximations to the tanh-curve, 
the exponential difference scheme itself is equivalent to first-order upwinding for 
convection with physical diffusion switched off, for physical grid Peclet numbers 
larger than about 6. This, of course, is entirely appropriate for the model physical 
problem on which Equation (19) is based; for P A > 6, the solution is essentially a 
constant (<J> = 0) at all computed points, since the effect of the downstream boundary 
condition does not diffuse upstream strongly enough to significantly affect the 
nearest computed point (e.g., e~ 6 ~ 0.25%). But this (steady, source-free, one-dimen- 
sional, convection-diffusion) model problem is clearly very restrictive. For this 
reason, it is dangerous to use Equation (19), or various equivalent formulations, as a 
basis for a general-purpose computational-fluid-dynamics algorithm. This is clearly 
demonstrated in the following section. 

Source Terms, Multidimensionality and Transients 

To the extent that the fluid-dynamic problem being modelled differs from 
Equation (1), Hybrid and PLDS (or EDS) introduce inaccuracies due to the inherent 
artificial diffusion - for P. > 2 in the case of Hybrid, and for all P. values for 
PLDS. For example, the introduction of a source term into Equation (1), represent- 
ing either a real source term in one dimension or the effects of transverse transport 
in a multidimensional flow, results in a serious degradation in accuracy when these 
methods are used. Figure 3 shows a steady-state two-point boundary-value problem 
corresponding to the model equation 


u 


3(f) 

dx 


d% 

= D — + Six) 

dx 


using a piece-wise linear source term given by [7] 


( 22 ) 


Six) ---- 10 - 50 iAx 

for 

(Ax < 0.3 | 


S(x) -- 50 (Ax — 20 

for 

0.3 < tAx < 0.4 J 

(23) 

Six) := 0 

for 

(Ax >0.4 ) 


where i is a grid index measured from the left-hand end and x = iAx . 

, Solutions 


are shown using both Hybrid and PLDS for P A = 1, 2, 6, and 10. Note that, for P A = 
1 and 2, Hybrid (operating as second-order central differencing) gives adequate solu- 
tions; but PLDS is already clearly in error even at these small P A values. Clearly, 
for IL > 2, the Hybrid solution does not change, since the effective P/ is frozen at 
2. With PLDS, the corresponding "saturated” solution occurs at about P A = 6, corre- 
sponding to Figure 2. Note that this is due to the shape of the hyperbolic-tangent 
function and has nothing to do with any errors involved in the power-law approxima- 
tion of tanh. 


To demonstrate the effects of multidimensionality, Figure 4 shows the well- 
known oblique-step test, corresponding to steady two-dimensional convection and 
diffusion governed by the model equation 


7 


cH|> ckj) 

a — + u — 
ax dy 



a 2 ^ 


+ 


a* 


<*V\ 

-.2 / 
<*y 7 


( 24 ) 


with a constant converting velocity, v = ( u,v ), at an angle 0 to the (uniform square) 
grid, and a prescribed unit-step jump in <]> on the upstream boundary. As 0 is 
varied, the position of the boundary jump is adjusted so that the centerline of the 
transition passes through the center-point of the grid. Figure 5 shows analytical 
solutions [32] for P A (= |v|Ax/£>) = 100, at 0 = 30°, 45°, and 60°. The analytical 
solution omits streamwise diffusion and is therefore valid only at moderate-to-high 
va l ues [33]. At this grid Peclet number, both Hybrid and PLDS are operating as 
pure first-order upwinding for convection with physical diffusion omitted. As seen in 
Figure 6, this results in gross artificial diffusion in a direction transverse to the flow 
— a well-known phenomenon [4], Note the large absolute error, given by 


£ — ^ | 4 > — 4 > ,| 

- — ’ comp ^exact 1 (25) 

summed over all computed grid points on a 25 X 25 square mesh. The effects of 
varying 0 (by 1° increments) are summarized in Figure 7. Note the symmetry with 
respect to 0 = 45°, where there is a slight drop in error due to alignment across 
diagonal nodes; other undulations are due to similar interactions between the 
almost-discontinuous exact profile and the discrete nodal points, as the angle is 
varied. Clearly, except when the direction of the convecting velocity is closely 
aligned with the grid (0 -> 0°, 90°, . . .), the artificial diffusion of these methods at 
practical (i.e., large) grid Peclet numbers results once again in highly inaccurate 
predictions. 

finally, Figure 8 shows results of simulating the transient pure-convection 
equation (for u = const > 0) 


# ckf> 

* + “ to = ° (26) 

for an isolated sine-squared profile 20Ax wide (representing a "smooth” function) 
and a unit step, using Hybrid or PLDS (or EDS) — equivalent to pure first-order 
upwinding for this infinite-P A case - and an explicit forward time step, i.e., 


*- +l = 4>;‘- c (<»"- <_,) (27) 

At the time shown, the exact profiles have moved 45 mesh-widths to the right. 
Starting from exact initial conditions, the computed results again demonstrate the 
gross artificial diffusion corresponding, in the transient case, to a numerical diffu- 
sion coefficient of the form [30] 


D — uAx (1 — c) / 2 (28) 

nuni 

In the case shown, c = 0.45; from Equation (28), numerical diffusion is clearly worse 
at smaller Courant numbers. 
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w b ^ 0Ve slm P Ie ,. test c r ase s show, Hybrid and PLDS are essentially equiva- 

lent to first-order upwinding of convection with physical diffusion ignored for Zid 
Peclet numbers larger than 2 or about 6, respectively. This must cast serious do^bt 
Hirn^ 6 ^ of simulations using these methods for complex high-speed multh 

dimensional flows using practical (i.e., relatively coarse) meshes In addS it 
seems to make no sense to use sophisticated (and expensive) multiple equation 
turbulence models with these convection schemes - as is usually done - because 

exceeds^ (H^briS C T?° n f ^ ^ id Peclet (or Reynolds) number 

tWelv) omittpd frnm fh (PLDS) ’ l + he turb “ lent , transport terms are (actually or effec 
t nrhni TI U d fro ™ the momentum and scalar transport equations (and from the 
nnmoif t i t ^fr sp ° rt Rations themselves) - having been replaced by the artificial 
numerical diffusion inherent in the truncation error of the modelled convection 

mnSflV Ap > rt ,?? m f he questionable logic involved, one wonders how turbulence 
modelling itself has been affected by such convection methods - since ^ except as an 
expensive diagnostic to switch itself off, turbulent transport is LrgelyTZred In 

r„1at t tLdZ y d w" y a ) ,n <eXCePt near S °' id b0Undaries whe - i wan-?unc 

PROBLEMS WITH HIGHER-ORDER METHODS 

Spurious Oscillation of Symmetrical Schemes 

th j art . 1 1 fi j ial ^^sion associated with Hybrid and PLDS is widelv 

poplar tvf34rThZ'l d ° CUme . nt l d ' ‘L 36 met . hods “"«"■« to enjoy widespread 
y . X *■ This appears to be because simple alternatives with formallv 
higher-order accuracy, introduce other anomalies - usually in tlie forZofoverTZ, 

hnea^ 0 s°tebnitt eV Z e e 0SC ' llatl0nS “, whic . h m “- v lead to convergence problems or non- 
•^f a ,J fui-il ty ’ rh P™ mar y motivation for Hybrid and PLDS seems to be asso 
ciated with the suppression of deeply penetrating oscillations whfch occur wfth 
classical” second-order central differencing fnr P ^ 9 ‘ oT ? ■ occur with 
Equation HA) Pim,r 0 qiT\ “l, 7l rencin £ tor La > 2 > as discussed in relation to 
} 4 : ^ e ,- a sbows the one-dimensional source-term boundary- value 
problem, considered earher, at P A = 10 solved using second-order central different 

<=inhiti'n S ^f en ’ ^ be oscillations penetrate from the downstream boundary well into the 
solution domam. In fact, the penetration distance is directly proportSto p" m 

By contrast, using second-order upwinding for convection T301 tno-pthpr with 

differenCing f0r di,fuSi0n - writteZ"omTo/-™!umf(CV^ 


al i | i + . + J 

D 

~ A? 


+ 1 


24>, + 


2 + 4> ( _,)- l - - 2 4>, _ j + 4> l 2 ) 




- S(x) 


(29) 

problem^ as st^n R^rfsibKo? P° n ^0 ' resu < ts f» r / hia . ona - d i m ensional 
accuracy than Hvhrid rvrPT fiq s method also gives much better 

severe undershoot and overshoot problems develop when using second order unwind ’ 
ingin t u>o-dimensional computations, as will be shown s «“ nd ->>rder upwmd- 
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Figure 9(c) shows the same one-dimensional problem solved by the QUICK 
method [37], or third-order upwinding, which, being a control-volume formulation, 
takes the form (for u > 0) 


U 

Ax 


( 4 , ; + 1 + 4 > , ) - " <4> i+ , - 2 4>i + 




D 

Ax 2 


^i + 1 i> 


= S(x) 


(30) 


Typically, a few oscillations are generated near the downstream boundary jump. 
But the penetration distance does not continue to grow with P a ; instead, it saturates 
at about 5Ax as P.-+ » [37]. Note, the much better (third-order) accuracy in the 
source-term region. "There has been some confusion in the literature concerning the 
third-order accuracy of QUICK for convection; this is explained in Appendix A. 


Many researchers appear to have an unsubstantiated aversion to unsymmet- 
rical higher-order (upwind) methods, and prefer to work instead with higher-order 
central differencing. For example, the fourth-order central discretization in control- 
volume form is given by 


U 

Ax 


7 ^ + 2 -^ + 1 - 4 > £ + 4 > l _ 1 ) 


(<b + <}> .) — 

" **-l 16 



1 

-! 
Ax 2 ' 


1 

- — (cj 
24 


-4> 


(^-^.i)- — (4» i+1 - 3 * i +3<l> j _ 1 - 


= Six) (31) 


But, as seen in Figure 9(d), although the smooth-region accuracy is good, the down- 
stream oscillation problem is actually worse than second-order! Once again, pene- 
tration distance is proportional to P A . This is a characteristic of all higher-order 
central-difference schemes; the proportionality constant actually increases with the 
formal order of accuracy [38]. From this simple test, it is seen why central-differ- 
encing (of any order) is difficult to work with. At this point, second and third ( an d 
possibly higher) order upwind methods seem to offer better alternatives in terms of 
accuracy and stability. 


Nonmonotonicity of Higher -Order Multidimensional Upwinding 

To demonstrate multidimensional effects, Figures 10 and 11 show second-order 
upwind and QUICK-2D [13] solutions, respectively, of the oblique-step test for P. - 
100, to be compared with the analytical and first-order solutions of Figures 5 and 6. 
Figure 10 clearly demonstrates the large undershoot and overshoot problems of 
second-order upwinding; but note the decrease in error from that of first-order 
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upwinding ( = Hybrid = PLDS). At 45°, QUICK gives qualitatively similar results 
to those of second-order upwinding, but with slightly less undershoot and overshoot 
and steeper resolution corresponding to the higher-order accuracy, as reflected by the 
smaller absolute error. Note the inverse relationship of undershoot and overshoot 
between these two methods as 0 is varied. Figure 12 shows the angular variation 
(in 1° increments) of absolute error (at P A = 100) for both methods; cumulative 
contributions less than zero (undershoots) or greater than one (overshoots) are also 
shown; note the asymmetry between overshoots and undershoots - i.e., an overshoot 
at 0 becomes an undershoot at (90° - 0) and vice versa, in each case. Again, note the 
drop in error at 45° and the undular behaviour at other angles. The two-dimensional 
oblique-step-test results of this and the previous section have been obtained by 
explicit time-marching to steady state; the conservative CV explicit flux formulation 
is summarized in Appendix B. 

This simple test problem clearly demonstrates the nonmonotonicity problems 
associated with higher-order methods in their "unlimited” form. What is needed is a 
simple limiting strategy which maintains the good resolution of higher-order up- 
winding while eliminating undershoots and overshoots without introducing artifi- 
cial diffusion or destroying conservation. The next section shows how this can be 
achieved in a straight-forward manner using an extremely simple technique which 
can be extended to arbitrarily high-order accuracy. 

THE ULTRA-SHARP ALTERNATIVE 

Normalized Variable Diagram for the Universal Limiter 

Figure 13 shows the local behaviour of the convected variable, §{x,y,z), in the 
vicinity of a control-volume (CV) face, in a direction normal to that face. Depending 
on the direction of the convecting velocity (here shown to the right), label the three 
indicated node values: <J> D (downwind), (upwind), and (centrally located 

between the other two), as shown. In part (a) of the figure, local behaviour is mono- 
tonic. One necessary condition of the universal limiter is sketched in the figure, 
shown by the cross-hatched limits: the convected face-value, 4>p should lie between 
adjacent node-values in locally monotonic regions. Note that this includes the limit- 
ing case, 4) c = ( = 4^)- But if = fyy, an additional condition is necessary, 

namely: $/- = = $u )* i n this case, as shown in part (b) of the figure. Treatment of 

locally nonmonotonic behaviour is sketched in parts (c) and (d) of Figure 13. Some 
flexibility is possible here, provided the interpolation is consistent, in terms of conti- 
nuity, with the previous conditions. 

It is more convenient to summarize the limiter constraints in terms of normal- 
ized variables. Figure 14(a) shows local behaviour near a CV face in a direction 
normal to the face. Note that node-value labelling again conforms to the convecting 
velocity direction - in this case, to the left. Now define the normalized variable, $ , 
anywhere in the vicinity of this CV face as 




4 >(x,y,z) - by 


(32) 


Note, in particular, that ^ = 0 and = 1. Figure 14(b) gives the same infor- 
mation as 14(a), but in terms of normalized variables. Now the universal limiter 
constraints can be portrayed in the Normalized Variable Diagram (NVD) - i.e., the 
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rv 

( £c». V P la ?l e - Thi , s is shown in Figure 15. Symbolically, the universal limiter con- 
straints on <|y can be written 


with 


^ < 1 


for 0 < <J> C < 1 


(33) 


$ f = 0 at <j> c = 0 


(34) 


and with continuous extensions beyond > 1 and $ 
monotonically rising condition: d<J> > 0. In Figure i5 

given by second-order central-differencing for <|> c > 1 : 


< 0 , maintaining a 
, these extensions are 


— 1 + 0.5 (<{> c — 1) for <j> c > 1 
and by second-order upwinding for negative 4> c : 


(35) 


= l .5? c for ^ < 0 (36) 

but other functional relationships are possible. The dotted lines show extensions of 
Equations (35) and (36) into the monotonic range, 0 < < 1; note that they both 

pass through the "second-order” point C , located at (0.5, 0.75). To avoid nonunique- 
ness near $ c -+ 0 + , the boundary OB has a steep but finite positive slope. This 
introduces an additional constraint 


<ty < const 4> c near 4> c 0 + (37) 

where const = 0(100), for example. 

At each stage of a pseudo- time-marching (Appendix B) or iterative solution, the 
universal limiter constraints are applied as follows: 

(i) For each CV face, note the direction of the (current) normal velocity com- 
ponent, thereby identifying 4>[/ - 4> c > anc * 4> 0 • 

(ii) Explicitly compute a tentative high-order multidimensional upwind- 
biassed face value, 

(iii) Compute the corresponding normalized face-value, and the normal- 
ized adjacent upwind node-value, 4 > c . 

(iv) If the point (§> c , $,) falls within the triangular region of Figure 15, simply 
proceed to the next CV face. 

(v) If not, <$> f is limited to the nearest appropriate constraint boundary at the 
given <J> C - value. 

(vi) The unnormalized face-value is then reconstructed by 


* f = + *u (38) 

(vii) The same procedure is used for each CV face. 
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Then in a time-marching algorithm, the (limited) convective fluxes and correspond- 
mg diffusive fluxes (computed using second-order central differencing) are now 
available for the explicit update step. Alternatively, iterative alternating-direction 
implicit solution can be implemented by introducing the downwind weighting factor 
as now described. ’ 

The Downwind Weighting Factor 

,t. *. n . st ? ad of limiting the face variable directly, the Downwind Weighting Factor 
(U Wh ) is introduced as an auxiliary variable, thereby generating a compact implicit 
scheme suitable for tridiagonal solution methods. After explicitly computing the 
high-order multidimensional upwind-biassed estimate, 4 >^, define 


^ 4> c (39) 

Since this is the same as 


1 - $c (40) 

it is not difficult to see that, in terms of DWF, the universal limiter constraints 
corresponding to Figure 15, become 


0 < DWF <1 for 0 < ft < 1 
or 


(41) 


DWF < 


(const - 1 )ft 


a - <t> c ) 

in the monotonic region, with, for example, 


near 


4v 


(42) 


and 


DWF = 




2(1 - ft c ) 


for ^<0 


(43) 


DWF = 0.5 for ft^, > 1 (44) 

This is shown in Figure 16. Note that the point A , given by (1,1) in Figure 15, has 
been stretched out into a vertical line in Figure 16. Point C is now at (0.5, 0.5). 

Now rewrite the face value (in terms of the known DWF) as 


= DWF 4> l 80 + (1 - DWF) <t»£ BC (45) 
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where "TJ3C” stands for "to-be-computed” - in the next iteration of an implicit line- 
sweep update. Note that the convective flux at each CV face now involves only the 
adjacent upwind and downwind " TBC ” node-values; but the DWF in Equation (45) 
implicitly contains higher-order wide-stencil information (in addition to the univer- 
sal limiter constraints). Clearly, when DWF’s are computed for each CV face on a 
quadrilateral grid, the implicit update stencil consists of the central CV node-value 
itself, together with only the surrounding adjacent node-values - i.e., a compact 
(2N + l)-point stencil in N dimensions. This is the same compact stencil as used by a 
large number of 2D and 3D commercial and research codes of TEACH-like structure, 
based on tridiagonal line-sweep solution algorithms, thus immediately opening up 
the possibility of incorporating (in principle, arbitrarily) high-order nonoscillatory 
convection schemes into these well-established general-purpose elliptic solvers. 

Using the popular compass-point notation [2], the individual convected face 
values are computed as follows. Considering the west face, for example, the face- 
value obtained from Equation (45) will depend on the sign of the convecting velocity 
at the west face. If is the initial higher order estimate on the west-face, the 
DWF is first computed according to 


or 


DWK + = if u >0 

w <t>p - 4v 


(46) 


DWK 




4> W -4> P 


if u <0 

W 


(47) 


Then (the appropriate) DWF is limited according to Figure 16. The face-value used 
in the implicit update is then 


or 


<t> = DWF + 4>„ + (1 - DWF + )4> 

^ w tv r iv W 


for it — 0 
tv 


(48) 


d> = DWF' $ w + (1 -DWF - )4> d for u <0 (49) 

^ IV LU W tv r w 

In either case, 4> is a linear combination of its two adjacent (upwind and downwind) 
node values, witfi SGNCu^) as a parameter, 


= (50) 

Similarly for the other faces of this particular CV cell. This results in an update (in 
three dimensions) of the form 

a p$p = a ^W + a S^S + V^B + a E$E + a N$N + a T$T + b ( 51 ) 

which is, of course, identical to the form generated by first-order upwind or second- 
order central methods or combinations such as Hybrid, PLDS, or EDS [2]. Again, it 
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should be stressed that the higher-order multidimensional information and the non- 
oscillatory universal limiter constraints are implicitly contained in the DWF’s of 
each face rather than involving "outlying” node-values that are then lumped into the 
explicit source term. The limiter constraints inherent in the DWF’s guarantee non- 
oscillatory results, with stability and convergence properties similar to first-order 
methods (DWF = 0) - but without introducing artificial numerical diffusion. 

Nonoscillatory Multidimensional Results 

Figure 17 shows results of applying the universal limiter constraints to the 
QUICK-2D scheme for the oblique step test at P. = 100. These ULTRA-QUICK 
results should be compared to the unlimited QUICK results of Figure 11. Note that 
resolution remains sharp (reflecting the third-order accuracy) even though there are 
no overshoots or undershoots. Also note the concomitant reduction in error seen in 
Figure 17(d). But the great power of the universal limiter resides in the fact that it 
can be used with higher-order convection schemes. A dramatic increase in sharpness 
in the simulation of the near-discontinuity is seen in Figure 18, showing results for 
the same problem using ULTRA-fifth-upwinding for convection. As explained in 
Appendix C, this scheme uses (one-dimensional) fifth-order upwinding normal to 
each CV face, together with third-order (QUICK) transverse curvature terms; 
second-order-central terms are adequate for diffusion. As 0 is varied from 0° to 90°, 
the total error is never more than about 5.2. Considering the fact that the error is a 
cumulative value summed over 625 grid-points, these ULTRA-5th results are seen to 
be highly accurate. 

Cost-Effectiveness and Adaptive Stencil Expansion 

Now that it has been established that viable nonoscillatory higher order 
methods can be easily devised, a natural question to ask is: which is more cost- 
effective in terms of overall computer usage ( for a prescribed accuracy), a fine-grid 
computation using a low-order method, or a coarse-grid calculation using a higher- 
order scheme? Gaskell and Lau [27] compared first-order upwinding with the third- 
order QUICK scheme and a nonoscillatory version (similar to ULTRA-QUICK), 
using successive grid refinement on the two-dimensional oblique-step test. For a 
prescribed accuracy, first-order methods required such a fine grid, compared with 
the practical-grid third-order methods, that CPU time was larger by three orders of 
magnitude! This disparity would be even more impressive in three dimensions. 
Another obvious question is: what order is "optimal” in the above sense of cost- 
effectiveness? Clearly second-order upwinding (which involves the same stencil as 
QUICK or ULTRA-QUICK) is not. Higher order nonoscillatory methods such as 
ULTRA-5th require more operations (than ULTRA-QUICK) at each grid point; but 
perhaps this can be offset by the greater accuracy available on a coarser grid. 

Figure 19 gives some indication of optimality in the two-dimensional case. 
Part (a) of the figure shows global absolute error versus number of grid-points for 
several different methods applied to the infinite-P A oblique step- test for 0 = 45°. 
Part .(b) shows the corresponding CPU time (without cost-adjustment for storage). 
The results, corresponding to a prescribed accuracy from part (a), are cross-plotted in 
part (c) of the figure. Even on a logarithmic scale, the "cost” of the first-order method 
is off-scale. As expected, third-order methods are significantly better than second- 
order upwinding (because of better accuracy for essentially the same cost-per-grid- 
point); but fifth and higher-order methods are somewhat less cost-effective than 
third. On the other hand, for a given grid-size, Figure 18 shows that ULTRA-5th, in 
particular, can give dramatically sharper resolution of discontinuities than the 
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corresponding third-order scheme. This suggests a simple strategy of adaptive 
stencil expansion that is proving to be extremely cost-effective in both steady-state 
and time-accurate transient simulations [39]: 

Use third-order ULTRA-QUICK in "smooth” regions of the flow; 
then on the basis of some "non-smoothness” monitor, switch to 
a higher-order (ULTRA) scheme locally, as needed. 

Since smooth regions correspond to small values of the convected variables’ 
change in gradient (i.e., curvature), a suitable monitor at face (i + i), for example 
would be the absolute average "curvature” (normal to the CV face): 


CURVAV = 0.5|(4> i+2 -<j>. + 1 )-(<j>. -4,^)1 (52) 

This, of course, is also a good indicator of sudden changes in gradient near discon- 
tinuities. Thus, if CURVAV is below a prescribed threshold, the algorithm uses the 
optimally cost-effective ULTRA-QUICK scheme. Even when discontinuities are 
present, this will account for the bulk of the flow domain - since, by definition, dis- 
continuities (when tightly resolved) involve only a very few number of grid points in 
narrow isolated regions. Wherever CURVAV exceeds the threshold, the algorithm 
automatically expands to a wider stencil, such as ULTRA-5th, and uses this more 
accurate scheme in the appropriate local regions. In most cases of practical interest, 
the additional "cost” is imperceptible compared with a global ULTRA-QUICK 
calculation; however, resolution of near-discontinuities is dramatically enhanced. 
Additional thresholds can be used for further adaptive stencil expansion to invoke 
(in principle, arbitrarily) higher-order resolution, if desired. In addition, it is 
desirable to use a higher-order (ULTRA) scheme where the local absolute normal 
gradient” 


GRAD = |*. + 1 -*.| (53) 

across the CV face exceeds a given threshold. This strategy has been used up to 
ninth-order in inviscid transient calculations [39] in order to resolve extremely 
narrow pulses in the convected variable. For steady-state calculations, it appears 
that an adaptive ULTRA-QUICK/5th/7th-order convection scheme has a number of 
attractive attributes, including cost-effectiveness (high coarse-grid accuracy), 
reliability (excellent stability and convergence properties), and ease of coding The 
order-switching strategy is shown schematically in Figure 20. Results of the 
oblique-step test (at P A = 100, with second-order diffusion terms) are seen in Figure 
21 to be essentially exact for all practical purposes. Even higher-order adaptive 
stencil expansion could be developed for specialty applications. 

Variable Curvature Factor 

~ . Fo1 ! second- and third-order methods, the normalized convected CV face-value, 
%[> is a single-valued function of <J>£, and the curve representing this relationship in 
the NVD passes through the point C of Figure 15. Passing through C guarantees 
second-order accuracy; third-order accuracy (in a control-volume formulation) 
requires a slope of 3/4 at C (see Appendix A). In terms of unnormalized variables, 
such schemes can always be written 


4y = 0.5 (<t) 0 + 4> c ) - VCF • (<t> D - 2<t, c + 4> y ) 
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where VCF = VCF(4> c ) is the variable curvature factor. For third-order methods, 
transverse-curvature terms [13] should be added for consistency. Equation (54) can 
be written in terms of normalized variables, as follows 

% f = 0.5 (1 + $ c ) - VCF .(i - 2$ c ) (55) 

Rearranging gives 


VCF 


(1 +$ c -2$y) 
2(1 - 2$ c ) 


= VCF^J 


(56) 


when 4>/<l> c ) is known. Note that this expression is determinate at <J> C = 0.5 for 
curves passing through point C . It is thus a simple matter to convert the universal 
limiter constraints of Figure 15 into constraints on VCF; these are shown in Figure 
22. For reference, the figure also shows the third-order QUICK value, VCF = 1/8. 
ULTRA-QUICK is shown by the heavy solid lines. Second-order central-differencing 
corresponds to VCF = 0; whereas second-order upwinding has VCF = 1/2. Third- 
order accuracy for a given "nonlinear” scheme requires VCF = 1/8 at <J> C = 1/2. 


In an iterative procedure, VCF is computed from Equation (56) based on 
current values (of )• Then Equation (54) is written in terms of "to-be-computed” 
values 


4V 


n k iaJ bc , 

05 (<t>y + 4 > C ) 


upu , . TBC 0 .TBC . .TBC. 
-VUM^ -2<1> C + 4^ ) 


(57) 


implying a pentadiagonal matrix algorithm (PDMA) for line-sweep solution. This 
technique was originally developed by Gaskell and Lau [27] using a "curvature 
correction”, a , from the default value of 1/8; i.e., in terms of the present notation, 


a($ c ) = VCF($' c ) - 1/8 (58) 

Figure 23 shows another third-order scheme based on an exponential upwinding or 
linear extrapolation refinement of third-order upwinding known as EULER [28]. 
Note the alternative treatment in the nonmonotonic regions. 

For reference, a pentadiagonal compass-point algorithm proceeds as follows: 

(i) For each CV face, compute VCF + or VCF -, depending on SGNlu^). 

(ii) The VCF value may be under-relaxed to enhance convergence [27]. 

(iii) Then for a typical, say west, face 

K = (° 5 - VCF^ + (0.5 + 2VCF t ;)d v - VCF^4» WW for « w >0 (59) 

or 

~ (0 5 ~ VCF~)4> w + (0.5 + 2VCF“)4> p - VCF“ $ E for ^<0 (60) 

In either case, is a linear combination of four adjacent node-values in a direction 
normal to the face (two upwind and two downwind), with SGN(uJ as a parameter: 
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(61) 


♦„ = X •+»■♦*•+*: SGN <“„>! 

[Inclusion of transverse curvature introduces additional terms which can be handled 
explicitly.] Similarly for the other faces of this CV cell. This results in an overall 
three-dimensional update algorithm of the form 

a r$r = a ww*ww + a sA'S + “bb^bb + °ee^ee + a i nn^nn + 

+ a w <b w + a s + a B 4> B + a £ <t> £ + + cy J> r + b (62) 

which is solved most effectively by alternating-direction pentadiagonal line-sweep 
iteration [27]. Tridiagonal methods could, of course, be used, taking care to place 
outlying - and certain other [17] - terms in the explicit source term; this, however, is 
not recommended as it appears to be somewhat less felicitous than the straight- 
forward (and highly efficient) PDMA technique. 

Relationship to TVD Schemes 

All of the so-called "shock-capturing” and "total-variation-diminishing” (TVD) 
schemes currently used in many gasdynamic codes [24] can be immediately adapted 
to incompressible convective modelling. Since these methods use the same stencil as 
QUICK (without transverse curvature), they can be described within the framework 
of the variable-curvature-factor technique. They can also, of course, be represented 
as specific curves in the downwind-weighting-factor diagram or the normalized- 
variable diagram. For example, Figures 24 and 25 show three well-known schemes 
of this type in terms of the NVD and the DWF, respectively: Superbee (a supercom- 
pressive scheme), Minmod (a relatively diffusive scheme) [26], and the "harmonic” 
scheme devised by van Leer [40]. 

All shock-capturing and TVD schemes in common use conform to an overly 
restrictive limiter [25] which tends to make them more diffusive than necessary. In 
terms of normalized variables, this restriction is 


$ < 2$ c for 0 < $ c < 0.5 (63) 

Superbee, for example, can be described as a hybrid between second-order central 
and second-order upwinding, constrained by Inequality (63) - and by 1 in the 

remainder of the monotonic region. Figure 26 shows the NVD for a scheme similar 
in philosophy to Superbee but constrained only by the universal limiter; this might 
appropriately be called ULTRA-B. In fact, the basic idea of just such a method was 
originally proposed by Roe for gasdynamic methods [41] and, coincidentally, named 
"Ultrabee”. Figure 27 shows the 45° oblique-step test at P A = 100 using Minmod, 
the harmonic scheme, and ULTRA-B (or Ultrabee). Note the diffusive nature of 
Minmod; the harmonic scheme gives results similar to those of ULTRA-QUICK; by 
contrast, the supercompressive ULTRA-B scheme gives extremely sharp monotonic 
resolution - in fact, slightly sharper than the exact solution at this finite grid Peclet 
number. 
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Although the step-resolution properties of simple supercompressive schemes 
such as ULTRA-B are very impressive, it should be realized that the sharp resolu- 
tion is actually obtained by the local introduction of negative numerical diffusion, 
which is directly responsible for artificial steepening of discontinuities. Unfortun- 
ately, such methods have a tendency to artificially steepen all profiles - with con- 
comitant flattening of peaks (due to mass conservation). This draw-back is seen in 
Figure 28, which shows oblique convection without diffusion (at 9 = 45°) of a semi- 
elliptical upstream boundary profile simulated using ULTRA-B; the (infinite-P A ) 
exact solution is shown for reference. By contrast, Figure 29(a) shows the same test 
simulated by ULTRA-QUICK/5th/7th. Note that, in this case, the sudden change in 
gradient at the base of the profile is well captured without the extreme clipping 
produced by the supercompressive scheme. However, there is still some loss of 
resolution near the local maximum because of the action of the universal limiter. 
This problem, in turn, can be largely resolved by switching off the limiter in regions 
of true local extrema, as seen in Figure 29(b). The problem then becomes one of 
pattern-recognition: near physical local extrema, the limiter needs to be switched 
off; but it should remain active near discontinuous jumps or sudden changes in 
gradient to suppress spurious unphysical overshoots or undershoots. An automatic 
discriminator of this type has been constructed for one-dimensional steady-state [7] 
and transient [39] cases. In Figure 29(b), the limiter is simply artificially switched 
off for 4> ^ 0.5. Clearly, this strategy cannot be used in general. An automatic dis- 
criminator for multidimensional steady-state flow is currently under development. 

Rotating Velocity Field 

To get some idea of the effect of spatially varying velocity, pure convection 
(without diffusion) of boundary-specified profiles in a velocity field given by solid- 
body rotation represents an effective benchmark problem. In a two-dimensional 
conservative, finite-volume, formulation, it is convenient to have the stream- 
function values available at control-volume corners; then the CV face-average 
convecting velocity is simply the difference of the corner stream-function values 
across the face. The stream-function used in the proposed benchmark test is 


CO 

= f K* -x.) 2 + (y -;y c ) 2 l (64) 

where co o is a specified angular velocity and (x c , y c ) is the location of the rotation axis. 
The velocity field is shown schematically in Figure 30. 

Figure 31 gives the exact solution for a sharp nondiffusing (P A = °° ) comple- 
mentary-error-function step in boundary value at the lower left corner. Figure 32 
shows results for first-order upwinding (= Hybrid = PLDS = EDS in this infinite- 
P A case). Figures 33 through 35 show unlimited forms of second-, third-, and fifth- 
order upwinding. Figures 36 and 37 give Minmod and ULTRA-B results. Finally, 
Figure 38 shows the ULTRA-QUICK/5th/7th results. These figures and the corre- 
sponding global errors listed in the captions are self-explanatory. 
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CONCLUSION 


^ In terms of normalized variables, first-order upwinding is represented by <J>, = 
fc which, in the NVD, is a straight line passing through (0,0) and (1,1). It thus con- 
forms (marginally) to the universal limiter constraints (Figure 15) and, thereby, 
produces nonoscillatory results. As shown by Godunov [42], it is the only "linear” 
scheme (i.e., in the present terminology, $ f is a linear function of ) which 
possesses this property, as can easily be seen from Figure 15. _TVD schemes are 
based on single-valued nonlinear relationships between and <]> r but also require 
the corresponding NVD curves to pass through (0,0) and (1,1), with certain other (in 
some respects, unnecessary) restrictions. If the NVD curve passes through (0.5, 
0.75), the scheme is second-order accurate (or third-order if it passes through this 
point with a slope of 3/4); if not, the scheme isjmly first-order accurate. Schemes 
with generally larger values of <|y (at a given d> c value) tend to be more "compres- 
sive than schemes with lower values. Thus, Superbee is more compressive than 
ULTRA-QUICK, for example; whereas, Minmod is more diffusive. First-order 
upwinding, of course, is seen to be the most diffusive of all nonoscillatory schemes. 

Higher-order schemes cannot be represented by a single curve in the (IjL,, ^ r ) 
plane; however, they can still be constrained to pass through (0,0) and (1,1), which 
appears to be a basic condition for nonoscillatory results. In the monotonic range 
(0 5s <]L 5s 1), the downwind^constraint, 1 , guarantees uniqueness, whereas the 
upwind constraint, 4>, <J> C , eliminates stair-casing. JThe additional ad hoc 

constraint, <fy # = const , avoids nonuniqueness near $ c -» 0 The resulting 

triangular region in the NVD, Figure 15, constitutes the steady-flow universal 
limiter in the monotonic range. In the nonmonotonic regime, there is considerable 
flexibility, provided the constraints are consistent with passing through (0,0) and 
(1,1), and that d<j>^/d$ c remains positive and finite (for uniqueness). It is not 
necessary to specify a single-valued function in this regime; this is usually done 
purely as a matter of convenience. 

As has been shown, it is a simple matter to desig^i higher-order nonoscillatory 
convection schemes and to frame them in terms of simple explicit time-marching 
techniques or (by using the downwind weighting factor) traditional TDMA iterative 
algorithms. (The special class of second- and third-order schemes (including all 
currently used TVD schemes) can also be constructed using straight-forward PDMA 
iterative methods, using the variable curvature factor technique - essentially 
equivalent to Gaskell-and-Lau’s "curvature correction” method.] Existing codes 
based on a TEACH-like structure involving iterative line-sweep TDMA solvers can 
now be immediately upgraded to incorporate cost-effective higher-order nonoscilla- 
tory methods such as ULTRA-QUICK/5th/7th. Alternatively, because of efficiencies 
involved in vectorization and parallel processing, the simple explicit time-marching 
algorithm may be more attractive. In any case, there is certainly no longer any 
reason why grossly artificially diffusive methods, such as Hybrid, PLDS (or EDS), 
should be used under any circumstances. The availability of ULTRA-SHARP 
schemes, giving nonoscillatory high-accuracy results on coarse grids, means that 
computational fluid dynamicists will at last be able to focus their attention on 
physical modelling of complex fluid phenomena and on reliable practical three- 
dimensional simulations for analysis and design, without being plagued by spurious 
oscillations or nonconvergence or by the (less obvious but potentially more 
dangerous) debilitating inaccuracies attributable to massive artificial diffusion. 
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APPENDIX A 

Formal Accuracy of the QUICK Scheme 

of accuracy oftha mST''" 81 ” th<! 1 j‘ erature . M3] regarding the formal order 

£ o • y 1 th ( \ UICK scheme in terms of truncation error. The problem seems to 

wfth fh!f en throug ^. a conf usion of the point-simulation of the convective derivative 
^ h . th corresponding control-volume formulation. It is enough to considlr the 

fa 1 ,? v, SC i? f ane ‘ dim ensional constant- velocity flow using a & uniform grid At 
third (and higher) order, there is a difference between ^ 


d<j> 

(A.l) 

on the one hand, and the integrated control-volume (CV) conservative i 


: expression 




Ax 


(A. 2) 


vfl h r^hl^ he A S T bS ? riptS refer t0 the ri £ ht and left CV face values of the convected 
variable. A Taylor-senes expansion of (A.l) results in [30] 

/#\ _ j^i-Vi *i + r 3 *i +3 Vrti 

2Ax " +0(Ax3) ( a3) 

XE2SZS: f he ^rd-difference, assuming a positive converting 
y. nis could be split into two terms, resulting in a conservative CV form 


d<{> 

dx 




Ax 


- + 0(Ax j 


where 


(A. 4) 


^ - 2 ( 4> 1 + 1 + 4 >,» - l ~ ( 4 >, + I - 2 < t > 4 + 


(A. 5) 


and is obtained from <J> r * by decreasing all indexes by 1. 

gi ven^y^A^ ^ to^tWrd ^oH T does t /^ represent the control-volume expression 
fbout SeCV example , ma ^ 6 ^ values 
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. / Ax 

K — 4> + <j> ( 

i + 1 ^ r \ 2 

K*: 

(f )'*;»:(“)■* i*"(f)‘* 

(A.6) 


'*:*:( 


(A. 7) 

, / 3 Ax ' 

l 1 "( 

3Ax\ 2 1 .../3Ax\ 3 1 A <iv»/ 3A * V 


1 - 4> r - 4> r \ 2 j 


~2~ ) ~ 6 V 2 ) + 24 V 2 ) 

(A. 8) 


So that linear interpolation across the right face would give 


4 V= jWVn + ty = <^ + < ^(-^) + 0(A * 4) 


(A. 9) 


It is easily seen that the 1/6 factor in Equation (A. 5) does not completely cancel the 
<J> " term. In order to cancel this term completely, one needs the QUICK formulation 


4>? = J MUi + ty- 


(A. 10) 


giving 


(b^ = d> + — <{> Ax'* + 

16 


(A. 11) 


Thus, for the complete CV term, 

4> - 4> 




e 1 

+ 

Ax 16 


4> - <j> 




Ax 3 + 


A . (A.12) 

Ax Ax 16 \ Ax ) 

and, of course, in the limit Ax — * 0, the term in square brackets becomes <J) ( , so that 


4>? - ♦? 4> r - 4> e 3 

r - + 0(Ax 3 ) 


Ax 


Ax 


(A. 13) 


Since a conservative flux-based formulation should be consistent with the CV form, 
the QUICK simulation is the appropriate one to use. Note that several other 
discretization schemes can be represented in terms of a curvature factor 


+?= ~ (4^ + 1 + 4>P — CF • (4> t + j — 2<j) ( + 4> ( _ j) (A. 14) 

(with an analogous expression for negative convecting velocities). For ex .^ i y i P^j’ 
second-order upwinding has CF = 1/2, and Fromm’s method [44] has CF = 1/4 (and, 
of course, central differencing is equivalent to CF = 0). The distinction between 
Equations (A.5) and (A. 10) thus hinges on the factors 1/6 and 1/8, respectively, for a 
steady-state control-volume formulation, the 1/8-factor gives the appropriate third- 


22 



order accuracy. Interestingly enough, for the corresponding time-accurate control- 
volume formulation of unsteady convection-dominated flows, an additional factor of 
1/24 (together with a Courant-number-dependent term) must be added to the 1/8- 
factor, thereby reconstructing the value of 1/6. This was explained - but evidently 
not universally understood - over ten years ago, when the QUICK and QUICKEST 
algorithms were initially published [37]. 

In terms of normalized variables, Equation ( A.10) becomes 

= 0.5(1 - 4> c ) - 0.125(1 - 2$ c ) (A.15) 

or, more simply 

= 0.75 + 0.75 ($£ - 0.5) (A. 16) 

Similarly, Equation (A. 14) becomes 

$ D = 0.75 + (2CF + 0.5) ($”„- 0.5) (A.17) 

r 

For finite CF, any scheme of this form passes through point C of Figure 15, and is 
therefore at least second-order accurate, as seen from Equation (A.14). If CF = 1/8, 
then the slope S is given by 

S = 2CF 4- 0.5 = 0.75 (A.18) 

so that a steady-state control-volume scheme which passes through point C with a 
slope of 3/4 is third-order accurate, according to Equation (A.16). Table A.l sum- 
marizes several second-order schemes in comparison with the third-order QUICK 
algorithm. 


Table A.l 


Method 

CF 

S 

Order 

QUICK 

1/8 

3/4 

Third 

Eq. (A. 5) 

1/6 

5/6 

Second 

Fromm 

1/4 

1 

Second 

2nd-up 

1/2 

3/2 

Second 

Central 

0 

1/2 

Second 
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APPENDIX B 


Explicit Control- Volume Time-Marching Algorithm 


i • V*[ in £ compass-point (plus "top” and "bottom”) notation, fluxes are computed 
explicitly for the west, south and bottom faces. For example, for the west face 


F = c <p -a GRAD 

LU LV IV W L 

where the west-face Courant number is 


(B.l) 


c = u A/ I Ax 

UJ w 


and the corresponding diffusion parameter is 


a ~ D At I A* 

W w 


(B.2) 


(B.3) 


The face- val ue is first computed explicitly, using (in general) a wide upwind-biassed 
multidimensional stencil; then the universal-limiter constraints are applied - again 
explicitly - to giye 6 . The gradient term is appropriately taken to be a second-order 
approximation simulating the diffusive flux; i.e., 


GRAD = I — ) Ax 
" ' & /«, 


- < 1 > 


w 


(B.4) 

[A higher order treatment of diffusion is not necessary because, for sharp changes in 
^ rac * 1 ^P 1 ^’ “fusion is very small, whereas when diffusion is large, gradients vary 
smoothly ] Then the update algorithm consists of a FORTRAN over-write assign- 
ment statement inside a multiple DO-loop: 

Set: 


k) 


k) + FJi.j, k) + k) + F b (i,j, k) 


- FJi+ 1 ,j, k ) - FUJ+ 1 , *) _ F b (i,j, k + l) + At SUJ, k) (B.5) 

where conservation has been used — e.g., for the east face at (i f j, k), 

F{i,j,k) = F w (i+\,j,k) (B.6) 

and so on Because of nonlinearities inherent in the limiter, it is not possible to give 
precise stability restrictions on the time-step. However, experience suggests that 
local Courant numbers near 0.3 appear to give satisfactory results under high- 
convection conditions. 

r ™e steady-state limiter portrayed in Figure 15 needs to be slightly modified 
lor the time-marching algorithm. Rather than using the ad hoc restriction given by 

accurate S limiter S [*29] ^ ^ & m ° re precise condition corresponding to the time- 
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rv 

s 4> c /c 


for 0 < <j> < 1 


(B.7) 


(in addition to requiring $ ^ 1 in the same range) where c is the appropriate 
normal-direction Courant number at the CV face in question. Without this con- 
straint, (very) slight undamped temporal oscillations may occur at isolated points 

neverteT e aoh S e C d° ntlnUlt ' eS a ‘ ^ P ‘ value * at such P» lnts ' a ^ady state might 

.. A1 ^P u ^ h onl y r } r st-order accurate in time, Equation (B.5) can give some indi- 
Frnm^nraM tempo r a 1 evolution of the solution from the assumed initial conditions, 
a P r f c fi ca lP° in t of view, this can indicate the total time necessary for reaching 
a steady state. One simply observes how long it takes for initial transients near an 

march hd forT to^wTlV ' ° f 3 f ownstrea . , 1 n boundary; then continued time- 

marching tor 2 to 3 times this time-scale is usually enough to cause reflections and 

devised"" 51 *" 6 aWay ' ° f C ° UrSe> m ° re prec ‘ se sto P pi "8 criteria can easfly“e 


APPENDIX C 

Higher-Order Multidimensional Fluxes 

EDS d ^cond^rde; r 1 ? ethods “ [ ncludin g combinations such as Hybrid, 

face 5 ’ are^nHrJtnT TVD ~. the multidimensional fluxes at each CV 

2Sn5^J d - 1 ?? th r ir norma ]-direction one-dimensional counterparts For 
6r u P™ lndm ?> for example, the left-face convective flux is given by upwind 
linear extrapolation in the normal direction, ** 

CFLUXL (ij, k) = CXL-(1.5<|> i _ l j k -0.5^ ^ (C.l) 

J e t r jf a XL 1S a PP r opriate left-face normal (x-component) Courant number Note 
that, to a consistent order, no "transverse curvature” terms are included 

exam^e,*the*QUICK iefW^»^onvectWeflu^is 1 gWen^^ VerSe C “ rVatUre; th “ S ’ f ° r 


where 


CFLUX L (i j k) = CXL • <b^ 


Q 


(C.2) 


using the definition 


4> 


<t> UN - l CURVLN + — CURVTL 

LJIN g 24 


(C.3) 


2 +<p l -l.j.k ) (C.4) 

and where, if CXL ^ 0, the normal and transverse "curvature” terms are given by 


CURVNL = ~ 2 < t > l _ 1 7 t 


■-2 j,k’ 


(C.5) 


25 


and 


CURVTL = (♦ i _ 1J+1 , 4 -2* j _ 1J ,* + *i-u-i > * ) 

+ ^i-lj,k+l ~ + U.k-J 

respectively; whereas, if CXL < 0, 

CURVNL = (* l+lj ' k - 2 * i j'k + * >i 


(C.6) 


(C.7) 


and 

CURVTL = (4> u+lt *- 2, l , ij,* + ^-u-i 1 * ) + ( *iJ.*+i- 2 *U> + 4> w.*-‘ ) (C ‘ 8) 

For higher-order convective fluxes, other transverse terms could be included. 
However, numerical experimentation has shown that 

(i) whereas inclusion of the third-order CURVT terms provides a significant 
increase in accuracy, 

(ii) higher-order transverse terms can actually reduce accuracy in some cases and 
are therefore omitted, and 

(iii) the coefficients of higher-order normal terms can be optimized to give better 
accuracy than that obtained from values using formal interpolation formulas. 

Thus, for example, the "fifth-order” convected value (at the left face) used in 
this report is given by 


• ( 5 ) i 

= 4 >, 


- - CURVAV + — FOURTH + -J- CURVTL 
6 128 24 


‘'f “ '♦'UN 

where CURVAV is an "average normal curvature” across the left face 

CURVAV = + 

FOURTH is the upwind-biassed fourth-difference normal to the face 


(C.9) 


(C.10) 


(C.ll) 


FOURTH = («j>. 4 l j k - 4* iJtk + 6 $,_i j, k - 4<l> ,— 2 j.k + 

(for CXL > 0), and CURVTL is as defined previously. Note the use of 1/6 (rather 
than 1/8) in Equation (C.9)! 

The "seventh-order” formula used in this report is given by 

4> < / ) = 4> un - l CURVAV + ^ FRTHAV + ~ SIXTH + ^ CURVTL (C.12) 
where FRTHAV is the (symmetrical) average fourth-difference across the left face, 
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FRTHAV = (<fr. +3 - 3<J). + 2 + 2*. + 1 + 2<t>. - 3<t» (1 + <J>._ 2 ) (C.13) 

(suppressing/s and fe’s), and SIXTH is the upwind-biassed sixth-difference normal to 
the face 


SIXTH = ( 4>. +3 - 6<J>. +2 + 15<J> [ + 1 + 20 < 1 ). + 154>. , - 6<t». 2 + <J>._ 3 ) ( C . 14 ) 

Formal interpolation would indicate factors of 1/8 multiplying CURVAV and 5/1024 
multiplying SIXTH - together with additional transverse and cross-difference terms. 
The factors used in Equation (C.12) have been found to give slightly more accurate 
results. 

In all cases, diffusive fluxes are modelled using the simple second-order approx- 
imation for the normal gradient. For example, for the left face, 

/a*\ 

\dxlt Ax (C.15) 

As explained previously, higher order terms are not warranted since, when diffusion 
is large, profiles are smooth and Equation (C.15) is consistent with the appropriate 
(QUICK) treatment of convection [13]; whereas, when steep gradients occur - neces- 
sitating higher-order convective modelling and the use of the universal limiter — 
diffusion terms are very small, so that a simple diffusion model suffices. 
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(C) HYBRID: P^ = 6. 


(Q) PLDS: P A = 6. 



FIGURE 3. - STEADY-STATE BOUNDARY- VALUE PROBLEM WITH SOURCE TERM GIVEN BY EQUATION (23). 
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0)0 = 30°, £ = 50.6. 



cb> e = 45°, £ = 64.2. 



(C) e = 60°, <f = 50.6. 


FIGURE 6. - CORRESPONDING SOLUTIONS OF THE OBLIQUE-STEP PROBLEM 
USING FIRST-ORDER UPW I ND I NG— EQUIVALENT TO HYBRID, PLDS, OR EDS 
AT THIS PECLET NUMBER (P A = 100). $ IS DEFINED IN EQUATION (25). 



FIGURE 7. - VARIATION OF TOTAL ABSOLUTE ERROR, «?, WITH STREAM-TO- 
GRID ANGLE, 0, IN 1° INCREMENTS FOR PLDS AT P A = 100 (EQUIVALENT 
TO FIRST-ORDER UPWINDING EXCEPT NEAR 0° AND 90°). MAXIMUM £ IS 
65.9 AT MM 0 OR 46°. 
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FIGURE 8. - FIRST-ORDER-UPWIND SOLUTION OF THE INVISCID TRANSIENT 
CONVECTION PROBLEM FOR A UNIT STEP AND AN ISOLATED SINE-SQUARED 


PROFILE, AFTER 100 TIME-STEPS AT 



(a) SECOND-ORDER CENTRAL DIFFERENCING. 



(b) SECOND-ORDER UPWINDING. 

FIGURE 9. - STEADY SOURCE-TERM BOUNDARY 


C = 0.45. 



(d) FOURTH-ORDER CENTRAL DIFFERENCING. 


VALUE PROBLEM OF FIGURE 3 FOR P A = 10. 
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(C) 0 = 60°, & = 21.2. (C) 6 = 60°, $ - 13.8 

FIGURE 10. - RESULTS OF THE OBLIQUE-STEP PROBLEM FOR P A = 100 USING FIGURE 11. - RESULTS OF THE OBLIQUE-STEP PROBLEM FOR P A = 100 USING 
SECOND-ORDER UPWINDING. QUICK-2D. 









®u 



(a) ORIGINAL UNNORMALIZED VARIABLES. 



FIGURE 1A. - DEFINITION OF NORMALIZED VARIABLES ACCORDING TO EQUATION (32). 



FIGURE 15. - UNIVERSAL LIMITER CONSTRAINTS IN THE NORMALIZED VARIABLE 
DIAGRAM (NVD). POINT C IS LOCATED AT (0.5. 0.75). 
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PI SURE 16. - UNIVERSAL LIMITER CONSTRAINTS (CORRESPONDING TO FIGURE 15) IN TERMS 
OF THE DOWNWINDWEIGHTING FACTOR (DWF) . POINT C IS AT (0.5. 0.5). NOTE THE 
SINGULARITY AT$ C * 1. 
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(C) e = 60°/ <? = 9.2. 



<d>£ AS A FUNCTION OF 0 (MAXIMUM = 10.2 AT 43° OR 47°). 
FIGURE 17. - ULTRA-QUICK OBLIQUE-STEP RESULTS FOR P A = 100. 


(C) 0 = 60°, (?= 5.0 



0° 

45° 

C 


(d)<* AS A FUNCTION OF 0 (MAXIMUM = 5.2 AT 10° OR 80°). 


FIGURE 18. - ULTRA-5TH-UPWIND OBLIQUE-STEP RESULTS FOR P. 
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FIGURE 22. - UNIVERSAL LIMITER BOUNDARIES PORTRAYED IN TERMS OF THE VARIABLE CURVATURE 
FACTOR. NOTE THE SINGULARITY AT <D C = 0.5. THE STRAIGHT LINES THROUGH 0 AND A CORRE- 
SPOND TO THOSE IN FIGURE 15. THE DASHED LINE SHOWS THE QUICK VALUE: VCF = 1/8. 
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FIGURE 24. - NORMALIZED VARIABLE DIAGRAM FOR SUPERBEE 
(HEAVY PIECEWISE LINEAR GRAPH), THE HARMONIC SCHEME 
(DOTTED PARABOLIC CURVE), AND MINMOD (DASHED LINES). 

ALL SCHEMES FOLLOW Of = O c IN THE NONMONOTONIC REGIONS. 



FIGURE 25. - DOWNWIND WEIGHTING FACTOR VARIATION FOR THE 
SUPERBEE (HEAVY LINES), HARMONIC (DOTTED), AND MINMOD 
(DASHES) SCHEMES. 



FIGURE 28. - NORMALIZED VARIABLE DIAGRAM FOR THE ULTRA-B 
SCHEME. COMPARE WITH SUPERBEE IN FIGURE 24, NOTE THE 
SIGNIFICANT DIFFERENCE NEAR O c - 0 + . 
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FIGURE 29. - OBLIQUE CONVECTION OF THE ELLIPTIC PROFILE. 


FIGURE 30. - SCHEMATIC DIAGRAM OF THE ROTATING- 
VELOCITY-FIELD BOUNDARY-STEP PROBLEM. 



FIGURE 31. - EXACT SOLUTION OF THE ROTATING -VELOCITY PROBLEM FOR 
P A = oo AND A NARROW COMPLEMENTARY-ERROR-FUNCTION PROFILE JUMP 
AT THE UPSTREAM BOUNDRY CORNER. 



FIGURE 32. - FIRST-ORDER UPWIND ( = HYBRID =PLDS= EDS) RESULTS 
FOR THE ROTATING PROBLEM; <£ = 65.3. 
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FIGURE 33. - SECOND-ORDER -UPWIND RESULTS FOR THE ROTATING STEP 
PROBLEM: $ = 27.3. 



FIGURE 39. - UNLIMITED QUICK-2D RESULTS FOR THE ROTATING STEP 
PROBLEM; £ * 18.9. 



FIGURE 35. - UNLIMITED FIFTH-ORDER UPWIND RESULTS FOR THE ROTATING 
STEP PROBLEM; £ = 16.1. 



FIGURE 36. - MINHOD RESULTS FOR THE ROTATING STEP PROBLEM; & » 30.6. 
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FIGURE 37. - ULTRA-B RESULTS FOR THE ROTATING STEP PROBLEM; 


<?= 7.7. 



FIGURE 38. - ULTRA-3rd/5th/7th RESULTS FOR THE ROTATING STEP 
PROBLEM; $ = 9.1. 
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advance for numerical convective heat transfer, for example. But, as this report shows, second-order TVD 
schemes form only a small, overly restrictive, subclass of a much more universal, and extremely simple, 
nonoscillatory flux-limiting strategy which can be applied to convection schemes of arbitrarily high-order 
accuracy, while requiring only a simple tridiagonal ADI line-solver, as used in the majority of general-purpose 
iterative codes for incompressible flow and numerical heat transfer. The new universal limiter and associated 
solution procedures form the so-called ULTRA-SHARP alternative for high-resolution nonoscillatory 
multidimensional steady-state high-speed convective modelling. 
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